Tree canopy
Population-based Scenario: AI: Increase by 10% in all zip codes
Targeted
Scenario AII1: Increase by 10% in zip codes in the lowest 1/5th of current TC cover (i.e. <=20th pctile)
Scenario AII2: Increase by 10% in zip codes in the highest 1/5th of the Social Vulnerability Index (i.e. >80th pctile)
Scenario AII3: Increase by 10% in zip codes in the highest 1/5th of hospitalization burden (i.e. >80th pctile)
Proportionate-universalism
Scenario AIII1: Increase by 10% for bottom 1/5th of current TC cover… down to 2% for top 1/5th
Scenario AIII2: Increase by 10% for top 1/5th of SVI … down to 2% for bottom 1/5th
Scenario AIII3: Increase by 10% for top 1/5th of hospitalization burden … down to 2% for bottom 1/5th
Impervious surface cover
Population-based: Scenario BI: Decrease by 10% in all zip codes
Targeted
Scenario BII1: Decrease by 10% in zip codes in the highest 1/5th of current imperv cover (i.e. >80th pctile)
Scenario BII2: Decrease by 10% in zip codes in the highest 1/5th of the Social Vulnerability Index (i.e. >80th pctile)
Scenario BII3: Decrease by 10% in zip codes in the highest 1/5th of hospitalization burden (i.e. >80th pctile)
Proportionate-universalism
Scenario BIII1: Decrease by 10% for top 1/5th of current imperv cover … down to 2% for bottom 1/5th
Scenario BIII2: Decrease by 10% for top 1/5th of SVI … down to 2% for bottom 1/5th
Scenario BIII3: Decrease by 10% for top 1/5th of hospitalization burden … down to 2% for bottom 1/5th
Number of zip codes with scenarios by scenario type
hosp_all_long %>%
filter(measure=="irr") %>%
group_by(scenario_intervention,scenario_type_7_abbrev) %>%
summarise(n_scenario=n()) %>%
ungroup()
## # A tibble: 14 × 3
## scenario_intervention scenario_type_7_abbrev n_scenario
## <chr> <chr> <int>
## 1 Imp PB 1260
## 2 Imp PU 1 1260
## 3 Imp PU 2 1260
## 4 Imp PU 3 1260
## 5 Imp T 1 1260
## 6 Imp T 2 1260
## 7 Imp T 3 1260
## 8 Tree PB 1260
## 9 Tree PU 1 1260
## 10 Tree PU 2 1260
## 11 Tree PU 3 1260
## 12 Tree T 1 1260
## 13 Tree T 2 1260
## 14 Tree T 3 1260
Implausible or unusual data Are there zip codes with ratios above 2? This implies the counterfactual scenario has twice as many cases as baseline
library(knitr)
n_by_scenario=hosp_all_long %>%
filter(measure=="irr") %>%
group_by(scenario_intervention,scenario_type_7_abbrev) %>%
summarise(n_scenario=n())
hosp_all_long %>%
filter(measure=="irr") %>%
filter(value>2) %>%
group_by(scenario_intervention,
scenario_type_3,scenario_sub_type, scenario_type_7_abbrev) %>%
summarise(n=n()) %>%
ungroup() %>%
left_join(n_by_scenario,by=c("scenario_intervention",
"scenario_type_7_abbrev")) %>%
mutate(prop=n/n_scenario) %>%
kable(digits=4)
| scenario_intervention | scenario_type_3 | scenario_sub_type | scenario_type_7_abbrev | n | n_scenario | prop |
|---|---|---|---|---|---|---|
| Imp | Population-based | NA | PB | 25 | 1260 | 0.0198 |
| Imp | Prop. Univ. | Exposure | PU 1 | 13 | 1260 | 0.0103 |
| Imp | Prop. Univ. | Hosp. burden | PU 3 | 10 | 1260 | 0.0079 |
| Imp | Prop. Univ. | SVI | PU 2 | 9 | 1260 | 0.0071 |
| Imp | Targeted | Exposure | T 1 | 7 | 1260 | 0.0056 |
| Imp | Targeted | SVI | T 2 | 3 | 1260 | 0.0024 |
| Tree | Population-based | NA | PB | 9 | 1260 | 0.0071 |
| Tree | Prop. Univ. | Exposure | PU 1 | 2 | 1260 | 0.0016 |
| Tree | Prop. Univ. | Hosp. burden | PU 3 | 4 | 1260 | 0.0032 |
| Tree | Prop. Univ. | SVI | PU 2 | 3 | 1260 | 0.0024 |
Ratios below zero? This implies the counterfactual scenario yields a negative number (negative predicted hospitalizations)
hosp_all_long %>%
filter(measure=="irr") %>%
filter(value<0) %>%
group_by(scenario_intervention,
scenario_type_3,scenario_sub_type, scenario_type_7_abbrev) %>%
summarise(n=n()) %>%
ungroup() %>%
left_join(n_by_scenario,by=c("scenario_intervention",
"scenario_type_7_abbrev")) %>%
mutate(prop=n/n_scenario)%>%
kable(digits=4)
| scenario_intervention | scenario_type_3 | scenario_sub_type | scenario_type_7_abbrev | n | n_scenario | prop |
|---|---|---|---|---|---|---|
| Imp | Population-based | NA | PB | 20 | 1260 | 0.0159 |
| Imp | Prop. Univ. | Exposure | PU 1 | 15 | 1260 | 0.0119 |
| Imp | Prop. Univ. | Hosp. burden | PU 3 | 12 | 1260 | 0.0095 |
| Imp | Prop. Univ. | SVI | PU 2 | 9 | 1260 | 0.0071 |
| Imp | Targeted | Exposure | T 1 | 5 | 1260 | 0.0040 |
| Imp | Targeted | SVI | T 2 | 1 | 1260 | 0.0008 |
| Tree | Population-based | NA | PB | 5 | 1260 | 0.0040 |
| Tree | Prop. Univ. | Exposure | PU 1 | 3 | 1260 | 0.0024 |
| Tree | Prop. Univ. | Hosp. burden | PU 3 | 3 | 1260 | 0.0024 |
| Tree | Prop. Univ. | SVI | PU 2 | 3 | 1260 | 0.0024 |
Facet by type of intervention (impervious surfaces vs tree canopy) and by scenario type - Population-based, Proportionate Universalism, Targeted
facet_histogram_fun=function(df){
df %>%
ggplot(aes(value))+
geom_histogram()+
facet_grid(
rows=vars(scenario_type_7_abbrev),
cols=vars(scenario_intervention)
)
}
#IRR
hosp_all_long %>%
filter(measure=="irr") %>%
#Remove ratios above 2
filter(value<2) %>%
#remove negative ratios
filter(value>0) %>%
facet_histogram_fun()+
xlab("IRR")
#IRD
hosp_all_long %>%
filter(measure=="ird") %>%
filter(value<1e-6) %>% #Exclude outliers
facet_histogram_fun()+
xlab("IRD")
#n cases
hosp_all_long %>%
filter(measure=="n_cases_diff") %>%
facet_histogram_fun()+
xlab("n_cases_diff")
Note I took the median of each measure over all zip codes.
hosp_all_long %>%
filter(measure=="irr") %>%
group_by(measure, scenario) %>%
summarise(irr_med=median(value,na.rm=TRUE)) %>%
ungroup() %>%
left_join(lookup_scenario, by ="scenario") %>%
ggplot(aes(y=irr_med,x=scenario_type_7_abbrev))+
geom_col()+
labs(x="Scenario", y="Median IRR over ZCTA")+
facet_grid(
#Move facet down to x-axis
#https://stackoverflow.com/questions/67519146/bar-plot-with-named-groups-on-x-axis-in-ggplot2
cols=vars(scenario_intervention),
scales="free_x",
space="free_x",
switch="x"
)+
theme(panel.spacing = unit(0, units = "cm"), # removes space between panels
strip.placement = "outside", # moves the states down
strip.background = element_rect(fill = "white")
)
## `summarise()` has grouped output by 'measure'. You can override using the
## `.groups` argument.
hosp_all_long %>%
filter(measure=="ird") %>%
group_by(measure, scenario) %>%
summarise(ird_med=median(value,na.rm=TRUE)) %>%
ungroup() %>%
left_join(lookup_scenario, by ="scenario") %>%
ggplot(aes(y=ird_med,x=scenario_type_7_abbrev))+
geom_col()+
labs(x="Scenario", y="Median IRD over ZCTA")+
facet_grid(
#Move facet down to x-axis
#https://stackoverflow.com/questions/67519146/bar-plot-with-named-groups-on-x-axis-in-ggplot2
cols=vars(scenario_intervention),
scales="free_x",
space="free_x",
switch="x"
)+
theme(panel.spacing = unit(0, units = "cm"), # removes space between panels
strip.placement = "outside", # moves the states down
strip.background = element_rect(fill = "white")
)
## `summarise()` has grouped output by 'measure'. You can override using the
## `.groups` argument.
hosp_all_long %>%
filter(measure=="n_cases_diff") %>%
group_by(measure, scenario) %>%
summarise(n_cases_diff_med=median(value,na.rm=TRUE)) %>%
ungroup() %>%
left_join(lookup_scenario, by ="scenario") %>%
ggplot(aes(y=n_cases_diff_med,x=scenario_type_7_abbrev))+
geom_col()+
labs(x="Scenario", y="Median diff. in N cases over ZCTA")+
facet_grid(
#Move facet down to x-axis
#https://stackoverflow.com/questions/67519146/bar-plot-with-named-groups-on-x-axis-in-ggplot2
cols=vars(scenario_intervention),
scales="free_x",
space="free_x",
switch="x"
)+
theme(panel.spacing = unit(0, units = "cm"), # removes space between panels
strip.placement = "outside", # moves the states down
strip.background = element_rect(fill = "white")
)
## `summarise()` has grouped output by 'measure'. You can override using the
## `.groups` argument.
hosp_all_long %>%
filter(measure=="n_cases_diff") %>%
group_by(measure, scenario) %>%
summarise(n_cases_diff_sum=sum(value,na.rm=TRUE)) %>%
ungroup() %>%
left_join(lookup_scenario, by ="scenario") %>%
ggplot(aes(y=n_cases_diff_sum,x=scenario_type_7_abbrev))+
geom_col()+
labs(x="Scenario", y="Sum diff. in N cases over ZCTA")+
facet_grid(
#Move facet down to x-axis
#https://stackoverflow.com/questions/67519146/bar-plot-with-named-groups-on-x-axis-in-ggplot2
cols=vars(scenario_intervention),
scales="free_x",
space="free_x",
switch="x"
)+
theme(panel.spacing = unit(0, units = "cm"), # removes space between panels
strip.placement = "outside", # moves the states down
strip.background = element_rect(fill = "white")
)
## `summarise()` has grouped output by 'measure'. You can override using the
## `.groups` argument.